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Abstract 

In the 60s and 70s super-massive central objects (from now onwards SMOs) were thought 
to be the main source of active galactic nuclei (AGNs) characteristics (luminosities of L « 
1O 12 L ). The release of gravitational binding energy by the accretion of material on to 
an SMO in the range of 10 7 - 10 9 M© has been suggested to be the primary powerhouse 
(Lynden-Bell 1969). That rather exotic idea in early time has become common sense nowa- 
days. Not only our own galaxy harbours a few million-solar mass black hole (Genzel 2001) 
but also many of other non-active galaxies show kinematic and gas-dynamic evidence of 
these objects (Magorrian et al. 1998) The concept of central super-massive stars (SMSs 
henceforth) (M > 5 x 10 4 M Q , where M. is the mass of the SMS) embedded in dense stel- 
lar systems was suggested as a possible explanation for high- energy emissions phenomena 
occurring in AGNs and quasars (Vilkoviski 1976, Hara 1978), such as X-ray emissions (Bah- 
call and Ostriker, 1975). SMSs and super-massive black holes (SMBHs) are two possibilities 
to explain the nature of SMOs, and SMSs may be an intermediate step towards the formation 
of SMBHs (Rees 1984). In this paper we give the equations that describe the dynamics of 
such a dense star-gas system which are the basis for the code that will be used in a prochain 
future to simulate this scenario. We also briefly draw the mathematical fundamentals of the 
code. 



1.1 The gaseous model 

To go from stationary to dynamical models we use a gaseous model of star clus- 
ters in its anisotropic version. It is based on the basic assumptions that: the system can be 
described by a one particle distribution function, the secular evolution is dominated by the 
cumulative effect of small angle deflections with small impact parameters (Fokker-Planck 
approximation, good for large iV-particle systems) and that the effect of the two-body relax- 
ation can be modelled by a local heat flux equation with an appropriately tailored conduc- 
tivity. 

The first assumption justifies a kinetic equation of the Boltzmann type with the inclusion of 
a collisional term of the Fokker-Planck (FP) type: 

®L +Vr ?L +v - r 2L +Vf) K +v - K=( 5 l\ (id 

dt ' dr ' dv r dvg v dv v ySt J ^ 
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In spherical symmetry polar coordinates r, 6, <j) are used and t denotes the time. The vector 
v = (v,-),i = r,6,<f) denotes the velocity in a local Cartesian coordinate system at the spatial 
point r,6,<f). The distribution function / is a function of, r, t, v r , v 2 , + only due to spherical 
symmetry. By multiplication of the Fokker-Planck equation (1.1) with various powers of 
the velocity components we get up to second order a set of moment equations which is 
equivalent to gas-dynamical equations coupled with Poisson's equation: A mass equation, a 
continuity equation, an Euler equation (force), radial and tangential energy equations. The 
system of equations is closed by a phenomenological heat flux equation for the flux of radial 
and tangential r.m.s. kinetic energy, both in radial direction. 

1.2 Interaction terms for the star component 

We now introduce the interaction terms to be added to right hand of the star com- 
ponent equations. 

7.2.7 Equation of continuity 

In the paper by Langbein et al. (1990) they derive the interaction terms to be added 
to the basic equations of the gaseous model. According to them, the star continuity equation 
is no longer 

■ + -2-5- (r 2 p*ii*) = 0, (1.2) 



di r 2 dx 



but 



where the right-hand term reflects the time variation of the star's density due to stars inter- 
actions (i.e. due to the calculation of the mean rate of gas production by stars collisions) and 
loss-cone (stars plunging onto the central object). 

If /(v re i) is the stellar distribution of relative velocities, then the mean rate of gas production 
by stellar collisions is 

P*fc{v K \) 



( s -£i) =-( 

V 6t / coU A, 



/(VrelWVl (1-4) 



COII J\V„\\ XTcoll f C0ll 

/(v re i) is a Schwarzschild-Boltzmann distribution, 

,, , 1 ( (Vrel.r-M*) 2 V rel,t\ 

/(Vrel)= 2^TV^- eXP ( 4*} la}) (L5) 

As regards f c , it is the relative fraction of mass liberated per stellar collision into the gaseous 
medium. Under certain assumptions given in the initial work of Spitzer & Saslaw (1966), 
we can calculate it as an average over all impact parameters resulting in r m ; n < 2r >t and as 
a function of the relative velocity at infinity of the two colliding stars, v re i. Langbein et al. 
(1990) approximate their result by 



/c(Vrel) = 

with <7 co ii = 100. So, we have that 



( 1 + q co n y/ (Tcoll/Vrel) V re l > <7 C0 U 

Vrel < CTcoll, 
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/c(Vrel) : 



0.01 (7 co ll = V re l 
CT co u > V re l, 



fcoii is defined as the mean time which has passed when the number of stars within a volume 
V = £ • v re i • Af is one, where v re i is the relative velocity at infinity of two colliding stars. 
Computed for an average distance of closest approach r„ 



this time is 



n*V(f co u) = 1 = «*Ev re ifcoii- 



And so, 



fcoll : 



with 



1 + - 



2Gm» 



'min^rel 



(1-6) 



(1.7) 



(1-8) 



ofe) = 2ctJ is the stellar velocity dispersion and E a collisional cross section with gravitational 
focusing. 

The first interaction term is * 



~5t 



~P*- 



-erff 



Ocoll 



-erff 



Ocoll 



VV6. 



(1-9) 



(1-10) 



coll tc °U 

which, for simplification, we re-call like this 

Sp *\ - n Y 

St /coll 

Since the evolution of the system that we are studying can be regarded as stationary, we 
introduce for each equation the "logarithmic variables" in order to study the evolution at 
long-term. In the other hand, if the system happens to have quick changes, we should use 
the "non-logarithmic" version of the equations. For this reason we will write at the end of 
each subsection the equation in terms of the logarithmic variables. 

In the case of the equation of continuity, we develop it and divide it by p* because we are 
looking for the logarithm of the stars density, d\np*/dt = (1/ p^dp^/dt. The result is: 



1 l-M* 

dt dr Or 



2w* 



1 

P* 



Sp* 
~5t 



1 

+ — 

coll P* 



Sp* 
~~5t 



(l.H) 



7.2.2 Momentum balance 

The second equation has the following star interaction terms: 



<9w* <9m* GM r 1 dp r Pr~Pt 
- + u r — — I- — — + — 1-2- 



dt 



dr 



p* dr 



p*r 



8u* 
~~5t 



(1.12) 



drag 



The interaction term is due to the decelerating force at which stars that move inside the gas 
are subject to. Explicitely, it is 



Su* 
~~St~ 



1 



— ^drag (w* Wg) 
drag P* 

* In the paper there are two typos in the equation, the correct signs and factors are given here. 



(1-13) 
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where we have introduced the following definition: 

2 

Xdrag = -C D P*PgCioU (1-14) 

with of ot = a 2 + of + (m* - u g ) 2 

To the end of the calculation of the logarithmic variable version of the equation, we mul- 
tiply by p*r/p r : 

p„r ( du % \ GM r dlnp r p t r 

-7T- + M* + P* + ^T, + 2 (1 ) = -Adrag— (1-15) 

p r \ at J rp r omr p r p r 

1.2.3 Radial energy equation 

As regards the last but one equation, the interaction terms are: 

-^- + -2TT( r U„p r ) + 2p r — + -— +~2^-( r F r) = (1-16) 

at r l or Or 5 f re iax r l or r 



where 



x j — -"urag^ r , i x 
St J drag V 6t / coll 

In order to determine e we introduce the ratio k of kinetic energy of the remaining mass after 
the encounter over its initial value (before the encounter); A: is a measure of the inelasticity 
of the collision: for k = 1 we have the minimal inelasticity, just the kinetic energy of the lib- 
erated mass fraction is dissipated, whereas if k < 1 a surplus amount of stellar kinetic energy 
is dissipated during the collision (tidal interactions and excitation of stellar oscillations). If 
we calculate the energy loss in the stellar system per unit volume as a function of k we obtain 

e = /7 1 [l-fe(l-/ c )]. (1.19) 

We divide by p r so that we get the logarithmic variable version of the equation. We make 
also the following substitution: 

F r = 3p r v r 

F t = 2p t v, (1.20) 

The resulting equation is 

d\np r , „ ,d\np r „(du* dv r \ 2( „ „ p,\ 
1 - + (M» + 3v r )^-^ + 3 — ^ + — ^ +- M* + 3v r -2v,— + 



* / drag V * / coll 



8pr\ 2 ( SPr \ 

-f— =-2XdragO- r , -FT =-X co llP*<Jr £■ (1-18) 



at Or \ or or ) r \ p r 

4 2 ~fr _ 1 fS Pr \ , 1 fSpr 



5 freiax Pr \ St J dmg p r \ St J ^ 
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1.2.4 Tangential energy equation 

To conclude the set of equations of the star component with the interaction terms, 
we have the following equation: 

-5T + -2-5-(r «*Pf) + 2 -— +~2T-( r F t) + — = (1-22) 

at r L or r 5 f re iax r or r 



^ = -2X dmg af, f% ) = -X coll p*a t 2 e. (1.24) 



where 



df / drag V df / coll 

We follow the same path like in the last case and so we get the following logarithmic 
variable equation: 

din p, dlnp, d 4 

— — + («» + 2v t ) — - + —("* + 2v,) + -(u* + 2v t ) - (1 .25) 

ot or or r 

I 2 !" 1 _ 1 (5p,\ + 1 ftp, 



5 freiax Pt \ St J dmg p, \ St J coU 

1.3 The gaseous component and the interaction terms 

In this section we give the set of equations corresponding to the gaseous component 
as for their right hand interaction terms. 

1.3.1 Equation of continuity 

For the SMS the equation of continuity looks as follows: 

where, for the mass conservation, we have that, obviously, 



Sp g \ ( 5p* 



5t /coll V * / coll 



(1.27) 



We follow the same procedure as for the star continuity equation to get the equation in 
terms of the logarithmic variables: 

dlnp g du g dlnp g 2u g _ 1 / Sp g \ 

^T + ^ +M ^ + V-^UJ coU (L28) 

The interaction term is in this case 

lfS Pg \_lf 6p.\__p. 

Pg 
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1.3.2 Momentum balance 

We modify equation number (2.9) of Langbein et al. (1990) in the following way: 

d( PgUg ) dp, dus 
dt Ug dt +Pg dt ' {L - M) 

we substitute this equality in their equation, divide by p g (u g is the variable in our code) and 
make use of the equation of continuity for the gas component. Thus, we get the following 
expression: 

du e du e GM r 1 dp r 4-n ( 5u e \ 

To get the interaction term we use the mass and momentum conservation: 
' S Pg\ . ( S P 



8( PgUg )\ + (SM\ =0 (133) 



6t /coll V St /coll 



We know that 

IT 



= 0, (1.34) 

coll 



5(p g u g )\ f Su g , , , - 

= w*p*A co n = p g l — — I +u g X C0 \\p^. (1.35) 



thus, 



S* Aon ' " S \St /coU 
Therefore, the resulting interaction term is 

Su g \ p* 

— J = — Xconiut-Ug) (1.36) 

6t / coll Pg 

In the case of the stellar system 

F=\(F r +F t )= 5 -p*v* (1.37) 
By analogy, we now introduce F m d in this way 

^=H= 5 -p g v g , (1.38) 
where v g is per gas particle. 

--is <L39) 

As means to write the equation in its "logarithmic variable version", we multiply the 
equation by p g r/p g , as we did for the corresponding momentum balance star equation and 
replace H by 5/2p g v g , 

p„r I du„ du„\ GM r dlnp„ 5 K ext 

— — - + iig—± } + p g + - -- p g rvg = 

Pg \ dt g dr ) rp g , s dlnr 2 c Fs 8 

r 

— p*X C0 \\(u*-u g ) (1-40) 
Pg 
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1.3.3 Radiation transfer 

We get the radiation transfer equations by re-writing the frequency-integrated mo- 
ment equations from Amaro-Seoane et al. (2003): We divide their first equation by J and 
multiply the second one by 2c/(5p g v g ), 

ldlnJ 5 3 5 3/ndd-l 

" + T7 i^:S v sPg) + -rPs v g — "s - ( 1 + /Edd) 



dt 2Jdr 81 k ' Jr' 8 4 cr 

1 0\np g _ Kabs 

c dt ~~T 



(B-J) (1.41) 



dlnvg | glngg | 2c 1 d{Jf E d&) | 2c 3/ E dd-l ^ 2« g 
<9? i9f 5 p^Vy dr 5 r/^Vg r 

-2^^=-cK ext p s (1.42) 

Where we have substituted H = 5p g v g /2 and /i aDS and n ex t are the absorption and extinction 
coefficients per unit mass 

_ PgMT) f \ n 

K abs — ~ j K ext — P?(, K abs + K scatUj (i-^J) 

t> 

A(T) is the cooling function, B the Planck function and K scatt the scattering coefficient per 
unit mass. We have made use of dM r /dr = 4n 2 p, fm& = K/J, and the Kirchhoff's law, 
B v = j v /n v (j v is the emission coefficient). 

1.3.4 Thermal energy conservation 

It is enlightening to construct an equation for the energy per volume unit e = (p r + 
2p t )/2 which, in the case of an isotropic gas (p r = p t ) is e = 3p/2. For this aim we take, for 
instance, equation (I1.17> and in the term 2p r du*/dr we include now a source for radiation 
pressure, 2(p r + p m A)du„/dr and we divide everything by e so that we get the logarithmic 
variables. The resulting equation is 

dine , „ Sine 2 / v \ f 8e\ 1 f 6e\ ft tn 

The interaction terms for this equation are 



■■X dmg (a 2 + cr 2 + (u^-u g ) 2 ) (1.45) 



drag 

5e\ 1 2 2 2 

T7 I = T^coiiP*((CT f + cr,) e + (u*-u g ) -£a coU ) (1.46) 

/ coll 1 

1.3.5 Mass conservation 

The mass conservation is guaranteed by 

3 dM r 

- = P* + P S (1-47) 



4vr dr 
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1.4 A mathematical view of the code 

Our model has seven dependent variables, 

(p,u,v r ,v,,p r ,p t ,M) = x(r,t), 



(1.48) 



So we have a set of non-linear, coupled differential equations plus an initial model, which 
very often is a Plummer's model. To solve our equations we discretise them on a logarith- 
mic radial mesh with typically 200 logarithmically equidistant mesh points, covering radial 
scales over eight orders of magnitude; e.g. from 100 pc down to 10~ 6 pc, which is enough to 
resolve the system down to the vicinity of a massive black hole's tidal disruption radius for 
stars. 

An implicit Newton-Raphson-Henyey iterative method is used to solve for the time evo- 
lution of our system. Let F(x) be a column-vector for the seven equations: 



F1O1...X7) \ 



F(X): 



F 1 (x 1 ...x 1 ) J 
The solution to the equations is x me ; therefore, 
F(x true ) = Q. 

Suppose that x (1) is a close value to the solution x lm ^ 
thus, 

F(x (1) ) = F(x true + Ax) = G(x), 



(1.49) 



(1.50) 



(1.51) 



(1.52) 



where G(x) is an "error function" that contains the difference between the exact value and 
the approximation. Since we have assumed that x (1) is a close value to x true , 



<9F 

F(x true + Ax) m F(x tme ) + — Ax = G(x) , 
ox 

therefore, since F(x true ) = 0, 



Ax = G(x) (§)"', 

dF\(xi . . .xj)/dxi 



(1.53) 



(1.54) 



where 



<9F 

dx 



\ 



(1.55) 



dF 7 (xi...x 7 )/dx 7 j 
is a 7 x 7 matrix. We iterate the process until we reach 

Ax (n) 



< £ - 10" 



(1.56) 



8 



P. Amaro-Seoane and R. Spurzem 

where n stands for the n-iteration done. This gives us the termination of the iteration. Our 
results compare well with other studies using direct solutions of the Fokker-Planck equation 
or Monte Carlo models (Lightman & Shapiro 1977, Marchant & Shapiro 1980). Note that 
the Monte Carlo approach has been recently revisited and improved by Freitag & Benz 
(2001). In contrast to the other models the gaseous model is much more versatile to include 
all kinds of important other physical effects, such as the dynamics of gas liberated in nuclei 
by stellar evolution and collisions and its interaction with the stellar component. 
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